/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#ifndef UHM_MESH_MESH_HXX
#define UHM_MESH_MESH_HXX

namespace uhm {
  typedef class Node_*      Node;
  typedef class Element_*   Element;
  typedef class Mesh_*      Mesh;
  typedef class Matrix_*    Matrix;
  typedef class Scheduler_* Scheduler;

  namespace interf {
    typedef class Sparse_*    Sparse;
    typedef class Mumps_*     Mumps;
    typedef class Pardiso_*   Pardiso;
    typedef class WSMP_*      WSMP;
    typedef class Petsc_*     Petsc;
  }

  //extern void mesh_new( Mesh &m );
  //extern void mesh_delete( Mesh &m );

  // internal inline functions
  bool   mesh_valid( Mesh m );

  // ----------------------------------------------------------------
  // ** Mesh class
  class Mesh_ : public Object_<int> {
  private:
    int id_element, locker, symmetry, method;
  protected:

    std::map < std::pair<int,int>, Node_ > nodes;
    std::map < int, Element_ > elements;

    Scheduler_ scheduler;

    void _init( int id, int id_element );
    void _random_matrix( int is_spd );

  public:
    Mesh_();
    Mesh_( int id );
    Mesh_( int id, int id_element );
    virtual ~Mesh_();

    virtual bool disp();
    virtual bool disp( int mode );
    virtual bool disp( FILE *stream );
    virtual bool disp( FILE *stream, int mode );

    bool import_file(char *full_path);
    bool export_graphviz_hier(char *full_path, int is_leaf2root);
    bool export_sparse_pattern(char *full_path, char *ss, int is_fill_in);
    bool export_connectivity(char *full_path, int n_rhs);
    bool export_matrix(char *full_path, int n_rhs);

    bool export_matrix(interf::Sparse sp, int n_rhs);
    bool export_matrix(interf::Sparse sp, int generation, int n_rhs);
    bool import_matrix(interf::Sparse sp, int assemble,
                       std::vector<double> &rhs, int ldb, int n_rhs);

    // external solver interface
    bool export_matrix(interf::Mumps mumps, int n_rhs, int is_sym);
    bool export_matrix(interf::Pardiso pardiso, int n_rhs, int is_sym);
    bool export_matrix(interf::WSMP wsmp, int n_rhs);
    bool export_matrix(interf::Petsc petsc, int n_rhs, int is_sym);

    bool operator<(const Mesh_ &b) const;

    void reset();

    void set_symmetry(int sym);
    int  get_symmetry();
    
    void set_decomposition_method(int method);
    int  get_decomposition_method();

    bool is_elements_separated();
    bool is_nodes_numbered();

    void backup( Mesh backup );
    void check_reuse( Mesh backup );

    Scheduler get_scheduler();

    int  get_n_nodes();
    int  get_n_elements();

    void remove_all_nodes();
    void remove_orphan_nodes();

    void remove_all_elements();
    void remove_lower_elements( int gen );
    void remove_non_leaves();
    void remove_leaves();
    void remove_non_boundary_nodes_from_leaves();

    // single physics
    Node add_node ( int id, int ndof );
    Node add_node ( int id, int ndof, int p );

    // multi physics
    Node add_node ( std::pair<int,int> id, int ndof );
    Node add_node ( std::pair<int,int> id, int ndof, int p );

    // advanced user
    Node add_node ( std::pair<int,int> id, int ndof, int p, int kind );

    Node add_node ( Node n );

    Node find_node( int id );
    Node find_node( std::pair<int, int> id );

    bool remove_node( int id );
    bool remove_node( std::pair<int,int> id );

    Element get_root();
    Element add_element();
    Element add_element( int gen );

    void    insert_element ( Element e );

    Element find_element   ( int id );
    Element insert_element ( int id );
    bool    remove_element ( int id );
    
    void adjust_element_numbering();

    Element refine_element  ( int id, int n_children );
    Element refine_element  ( int id, int is_binary, int n_children );
    Element unrefine_element( int id );

    // solving sequence
    void lock();
    void unlock();
    int  is_locked();

    void create_matrix_without_buffer( int datatype, int n_rhs );

    // i want to separate create buffer between leaf and element
    void create_leaf_matrix_buffer();
    void create_element_matrix_buffer();
    void create_element_matrix_buffer(int is_schur);
    
    // old one same as element buffer creation
    void create_matrix_buffer();
    void create_matrix_buffer( int is_schur );

    void free_matrix();
    void free_matrix_buffer();
 
    void random_matrix();
    void random_spd_matrix();
    void triangularize();

    void copy_in  ( Element elt,
		    int datatype, int m, int n,
		    int *nods, int side,
		    void *buffer );
    void copy_in  ( Element elt,
		    int datatype, int m, int n,
		    int nod_disp, int *nods, int side,
		    void *buffer );
    void copy_in  ( Element elt,
		    std::vector< std::pair<int,int> > &nods, 
		    int side,
		    linal::Flat_ A );

    void copy_out ( Element elt,
		    int datatype, int m, int n,
		    int *nods, int side,
		    void *buffer );
    void copy_out ( Element elt,
		    int datatype, int m, int n,
		    int nod_disp, int *nods, int side,
		    void *buffer );
    void copy_out ( Element elt,
		    std::vector< std::pair<int,int> > &nods, 
		    int side,
		    linal::Flat_ B );

    void         set_rhs();
    double       get_residual();
    double       get_lower_triangular_norm();
    unsigned int get_n_dof();
    unsigned int get_n_nonzero_factor(int is_sym);
    unsigned int get_n_nonzero();

    void   estimate_cost( int method, int datatype, int n_rhs,
                          double &flop_decompose, double &flop_solve,
                          double &buffer );
    void   get_matrix_stat( int method, int datatype,
                            int delta_dim,
                            std::vector< int > &x,
                            std::vector< int > &dist_uhm,
                            std::vector< double > &dist_flop,
                            std::vector< double > &dist_buffer );

    void export_graphviz_lu_nopiv(char *full_path, int bmn);


    // ** I am not going to complete all ooc part. 
    //    It is not necessary now. Maybe forever...
    
    void decompose_with_free    ( int method );
    void decompose_without_free ( int method );
    void solve_1                ( int method );
    void solve_2                ( int method );
    void solve                  ( int method );
    void check_1                ( int method );
    void check_2                ( int method );
    void check                  ( int method );

    void decompose_with_free    ();
    void decompose_without_free ();
    void solve_1                ();
    void solve_2                ();
    void solve                  ();
    void check_1                ();
    void check_2                ();
    void check                  ();

    // ---------------------
    void chol_with_free();
    void chol_without_free();
    void chol_with_ooc();       // not done
    // ---------------------
    void solve_chol_1();
    void solve_chol_2();
    void solve_chol();

    void solve_chol_1_ooc();   // not done
    void solve_chol_2_ooc();   // not done
    void solve_chol_ooc();     // not done

    void check_chol_1();
    void check_chol_2();
    void check_chol();

    void improve_chol();       // not done
    // ---------------------
    void chol_piv_with_free();
    void chol_piv_without_free();
    void chol_piv_with_ooc();
    // ---------------------
    void chol_incpiv_with_free();
    void chol_incpiv_without_free();
    void chol_incpiv_with_ooc();
    // ---------------------
    void solve_chol_piv_1();
    void solve_chol_piv_2();
    void solve_chol_piv();

    void solve_chol_piv_1_ooc();
    void solve_chol_piv_2_ooc();
    void solve_chol_piv_ooc();

    void check_chol_piv_1();
    void check_chol_piv_2();
    void check_chol_piv();

    void check_chol_piv_1_ooc();
    void check_chol_piv_2_ooc();
    void check_chol_piv_ooc();

    void improve_chol_piv();      // not done
    void improve_chol_piv_ooc();  // not done
    // ---------------------
    void ldl_nopiv_with_free();
    void ldl_nopiv_without_free();
    void ldl_nopiv_with_ooc();   // not done
    // ---------------------
    void solve_ldl_nopiv_1();
    void solve_ldl_nopiv_2();
    void solve_ldl_nopiv();

    void solve_ldl_nopiv_1_ooc(); // not done
    void solve_ldl_nopiv_2_ooc(); // not done
    void solve_ldl_nopiv_ooc();   // not done

    void check_ldl_nopiv_1();
    void check_ldl_nopiv_2();
    void check_ldl_nopiv();

    void improve_ldl_nopiv();     // not done

    // ---------------------
    void ldl_piv_with_free();
    void ldl_piv_without_free();
    void ldl_piv_with_ooc();
    // ---------------------
    void ldl_incpiv_with_free();
    void ldl_incpiv_without_free();
    void ldl_incpiv_with_ooc();
    // ---------------------
    void solve_ldl_piv_1();
    void solve_ldl_piv_2();
    void solve_ldl_piv();

    void solve_ldl_piv_1_ooc();
    void solve_ldl_piv_2_ooc();
    void solve_ldl_piv_ooc();

    void check_ldl_piv_1();
    void check_ldl_piv_2();
    void check_ldl_piv();

    void check_ldl_piv_1_ooc();
    void check_ldl_piv_2_ooc();
    void check_ldl_piv_ooc();

    void improve_ldl_piv();      // not done
    void improve_ldl_piv_ooc();  // not done
    // ---------------------
    void lu_nopiv_with_free();
    void lu_nopiv_without_free();
    void lu_nopiv_with_ooc();   // not done
    // ---------------------
    void solve_lu_nopiv_1();
    void solve_lu_nopiv_2();
    void solve_lu_nopiv();

    void solve_lu_nopiv_1_ooc(); // not done
    void solve_lu_nopiv_2_ooc(); // not done
    void solve_lu_nopiv_ooc();   // not done

    void check_lu_nopiv_1();
    void check_lu_nopiv_2();
    void check_lu_nopiv();

    void improve_lu_nopiv();     // not done
    // ---------------------
    void lu_piv_with_free();
    void lu_piv_without_free();
    void lu_piv_with_ooc();
    // ---------------------
    void lu_incpiv_with_free();
    void lu_incpiv_without_free();
    void lu_incpiv_with_ooc();
    // ---------------------
    void solve_lu_piv_1();
    void solve_lu_piv_2();
    void solve_lu_piv();

    void solve_lu_piv_1_ooc();
    void solve_lu_piv_2_ooc();
    void solve_lu_piv_ooc();

    void check_lu_piv_1();
    void check_lu_piv_2();
    void check_lu_piv();

    void check_lu_piv_1_ooc();
    void check_lu_piv_2_ooc();
    void check_lu_piv_ooc();

    void improve_lu_piv();      // not done
    void improve_lu_piv_ooc();  // not done
    // ---------------------
    void qr_with_free();
    void qr_without_free();
    void qr_with_ooc();       // not done
    // ---------------------
    void solve_qr_1();    // done
    void solve_qr_2();    // done
    void solve_qr();      // done :  not verified

    void check_qr_1();    // not done
    void check_qr_2();    // not done
    void check_qr();      // not done

    void improve_qr();    // not done
    // ---------------------
    friend bool mesh_valid( Mesh m );
    friend bool build_tree_var_1( Mesh m );
    friend bool build_tree_var_2( Mesh m, int nparts );
    friend bool build_tree_var_3( Mesh m );
    friend bool build_tree_var_4( Mesh m );
    friend bool build_tree_var_5( Mesh m );
    friend bool partition_mesh_var_1( Mesh m, int nparts );
    friend bool partition_mesh_var_2( Mesh m, int nparts );
    friend class Scheduler_;
  };
  // ----------------------------------------------------------------
  // ** Definition
  inline void Mesh_::_init( int id, int id_element ) {
    this->cookie     = UHM_MESH_COOKIE;
    this->id         = id;
    this->id_element = id_element;
    this->symmetry   = UHM_UNSYMMETRY;
    this->locker     = 0;
    this->method     = UHM_LU_PIV;
  }
  inline bool Mesh_::operator<(const Mesh_ &b) const { 
    return (this->id < b.id); 
  }
  inline bool mesh_valid( Mesh m ) { 
    return (m && m->cookie == UHM_MESH_COOKIE); 
  }

}

#endif
